#ifndef EXPDISK_H
#define EXPDISK_H

#include <stdio.h>

#include "MACRO.H"
#include "Vector.h"
#include "Particle.h"




typedef unsigned int uint;
typedef long unsigned int luint;

using namespace std;


class ExpDisk {
	public:
		

		double m;
		double r0;
		double z0;
		double rmax;
		double zmax;
	
		
		ExpDisk() :
		m(1),
		r0(2.5),
		z0(0.5),
		rmax(100),
		zmax(100)
		{}
		
		
		// ok
		double rho(double r, double z) {
			return m/(4*M_PI*z0*r0*r0)*(1-pow(tanh(z/z0),2))*exp(-r/r0);
		};
		
		
		template <typename T> 
		void set_rm(vector<T> &p) {
			
			printf("ExpDisk::set_pos()... "); fflush(stdout);
		
			srand48(222);

			double r, z, r_old;
			
 			for(int i = 0; i<p.size(); i++)  {
			
				do {
				
					double q = drand48();
					z = z0 / 2 * log(q / (1 - q));
					q = drand48();
					r = 1.0;
					
					do	{
					
						double f = (1 + r) * exp(-r) + q - 1;
						double f_ = -r * exp(-r);

						r_old = r;
						r = r - f / f_;
					
					} while(fabs(r - r_old) / r > 1e-7);

					r *= r0;

				} while (r > rmax || fabs(z) > zmax);	//countr++; else i++;
				
				double phi = drand48() * M_PI* 2;
				
				p[i].R[0] = r * cos(phi);
				p[i].R[1] = r * sin(phi);
				p[i].R[2] = z;
				
			}
		
			for(int i=0; i<p.size(); i++)
				p[i].m = m/p.size();
			
			printf("done\n");
			
		}
		
		void print() {
			
			printf("ExponentialDisk::print()\n");
			printf(" ├ m = %g\n", m);
			printf(" ├ r0 = %g\n", r0);
			printf(" ├ z0 = %g\n", z0);
			printf(" ├ rmax = %g\n", rmax);
			printf(" └ zmax = %g\n", zmax);
			//printf(" └ G = %g\n\n", G);
				
		}
		
};

#endif
